clear all;clc;close all;
% Последовательная RL цепь, выход - напряжение с сопротивления
L = 1e-4;
R = 50;

w_max = R/L * 2*pi* 10;
dw = w_max/1000;
w = 0:dw:w_max-dw;

a = [L/R 1];
b = 1;

% ПФ
% figure(1)
% freqs(b,a,w);

H = freqs(b,a,w);
% 
% figure(2)
% subplot(2,1,1)
% plot(w/2/pi,abs(H));
% xlabel('f,Hz');
% ylabel('|H(f)|');
% subplot(2,1,2)
% plot(w/2/pi,180/pi*unwrap(angle(H)));
% xlabel('f,Hz');
% ylabel('|arg(H(f))|,deg');

%ИХ
sys = tf(b,a);
[h,t] = impulse(sys);

% figure(3);
% plot(t,h);

% Дискретный аналог
% [z,p,k] = tf2zp(b,a);
% [dz,dp,dk]=bilinear(z,p,k);

Fs = 10*R/L;
[db,da]=bilinear(b,a,Fs);
[Hz,wz]=freqz(db,da);

% figure;
% plot(w/2/pi, abs(H),Fs*wz/2/pi, abs(Hz));
% legend('analog','digital')

% Пропускаем шум
dt = 1/Fs;
t_mod = 10000 * dt;
t = 0:dt:t_mod-dt;

stdn = 20.07;
n = stdn*randn(1,length(t));

y = filter(db,da,n);

% p=plot(t,[n;y]);
% set(p(2),'LineWidth',2,'Color','r');

s_n = fft(n)/length(n);
f = 0:1/t_mod:(1/dt - 1/t_mod);
s_y = fft(y)/length(y);

mn = sqrt(mean(abs(s_n).^2))*2.5;

s_y = abs(s_y)/mn;
s_n = abs(s_n)/mn;

plot(f,[s_n;s_y],w/2/pi,abs(H),'LineWidth',3);
xlim([0 Fs/2]);

